This class controls equilibration of atomic charges in the system.
It is possible for the user to define the charges of atoms in ASE. If a system exhibits charge transfer, polarization, charged defects etc., one may not know the charges beforehand or the charges may change dynamically during simulation. To handle such systems, it is possible to let the charges in the system develop dynamically.
Since charge dynamics are usually much faster than dynamics of the ions, it is usually reasonable to allow the charges to equilibrate between ionic steps. This does not conserve energy exactly, however, since the charge equilibration drives the system charge distribution towards a lower energy. The energy change in charge redistribution is lost unless it is fed back to the system.
Special care must be taken when setting up links between the atomic structure (ASE Atoms), the calculator (Pysic), and the charge relaxation algorithm (ChargeRelaxation). While some of the objects must know the others, in some cases the behavior of the simulator changes depending on whether or not they have access to the other objects.
The atoms and the calculator are linked as required in the ASE calculator interface: One can link the two by either the set_atoms() method of Pysic, or the set_calculator method of ASE Atoms. In either case, the atomic structure is given a link to the calculator, and a copy of the structure is stored in the calculator. This must be done in order to do any calculations on the system.
Also the relaxation algorithm has to know the Pysic calculator, since the relaxation is done according to the Potential interactions stored in the calculator. The algorithm can be made to know the calculator via the set_calculator() method of ChargeRelaxation. By default, this does not make the calculator know the relaxation algorithm, however. Only if the optional argument reciprocal=True is given the backwards link is also made. Pysic can be made to know the relaxation algorithm also by calling the method set_charge_relaxation(). Unlike the opposite case, by making the link from the calculator, the backwards link from the relaxation algorithm is always made automatically. In fact, even though linking an algorithm to a calculator does not automatically link the calculator to the algorithm, if a different calculator was linked to the algorithm, the link is automatically removed.
This slightly complicated behavior is summarized as follows: The algorithm should always have a link to a calculator, but a calculator need not have a link to an algorithm. If a calculator does link to an algorithm, the algorithm must link back to the same calculator. Clearly one does not always want to perform charge relaxation on the system and so it makes sense that the calculator need not have a link to a charge relaxation algorithm. If such a link does exist, then the relaxation is automatically invoked prior to each energy and force evaluation. This is necessary in simulations such as molecular dynamics (MD). A charge relaxation can be linked to a calculator in order to do charge equilibration, but if one does not wish to trigger the charge relaxation automatically, then it is enough to just not let the calculator know the relaxation algorithm.
The atomic structure cannot be given a link to the relaxation algorithm since the charge relaxation is not part of the ASE API and so the atoms object does not know how to interact with it. In essence, from the point of view of the structure, the charge relaxation is fully contained in the calculator.
The charge relaxation algorithm always acts on the structure contained in the calculator. The atomic charges of this structure are automatically updated during the relaxation. Since the calculator only stores a copy of the original structure, the original is not updated. This may be desired if, for instance, one wishes to revert back to the original charges. However, during strucural dynamics simulations such as MD, it is necessary that the relaxed charges are saved between structural steps. This is a problem, since structural dynamics are handled by ASE, and ASE invokes the calculation of forces with the original ASE Atoms object. Therefore, if the relaxed charges are not saved, charge relaxation is always started from the original charges, which may be very inefficient. In order to have also the original structure updated automatically, the charge relaxation can be made to know the original structure with set_atoms(). Note that the structure given to the algorithm is not used in the actual relaxation; the algorithm always works on the structure in the calculator, which may be different. The given structure is merely updated according to the calculation results.
Sometimes the most efficient way to optimize a system is to combine several algorithms or several parameterizations of a single algorithm in a sequence. First, one can run a robust but less efficient algorithm to find a rough solution and follow with a more efficient method to pinpoint the optimum at high precision. For a plain relaxation, this can be done manually. However, if the optimization needs to happen automatically, e.g., between molecular dynamics steps, one needs to combine several ChargeRelaxation objects together. This can be done through the method add_relaxation_pipe().
Adding a piped ChargeRelaxation object simply calls the charge_relaxation() methods of both objects one after another when relaxation is invoked. Several objects can be piped by recursive piping. For instance, the following:
rel1 = ChargeRelaxation(...)
rel2 = ChargeRelaxation(...)
rel3 = ChargeRelaxation(...)
rel1.add_relaxation_pipe(rel2)
rel2.add_relaxation_pipe(rel3)
creates a pipe rel1-rel2-rel3, where rel1 is executed first and rel3 last. Creating looped pipes where the same relaxation appears twice is not allowed.
By default, when a pipe is created, the calculator and structure objects of the first relaxation are copied also to all the piped ChargeRelaxation objects (including recursively defined ones).
Similarly to ASE molecular dynamics, also the ChargeRelaxation allows one to attach callback functions to the dynamics. That is, one can have a specific functions be automatically invoked at given intervals. This can be used for instance for printing statistics or saving data during the optimization run. For instance, the following example saves and plots the charge of the atom with index 0:
system = Atoms(...)
calc = Pysic(...)
rel = ChargeRelaxation(calculator = calc, system = system)
charge0 = []
def save_charge():
charge0.append(calc.get_atoms().get_charges()[0])
rel.add_observer(save_charge)
rel.charge_relaxation()
import matplotlib.pyplot as plt
plt.plot(np.array(range(len(charge0))), thecharges)
plt.show()
Below is a list of the charge relaxation methods currently implemented.
Assigning an inertia, \(M_q\), on the atomic charges, \(q_i\), we can describe the system with the Lagrangian
where \(m_i, \mathbf{r}_i\) are the mass and position of atom \(i\), respectively. The last term is a Lagrange multiplier corresponding to the constraint of fixed total charge, i.e., \(\sum_i q_i = Q_\mathrm{tot}\) being constant. The total potential energy \(U\) is a function of all charges and positions.
The equations of motion for this system are
In the charge equation, the Lagrange multiplier can be shown to equal the average electronegativity of the system, \(\nu = \bar{\chi}\), and the derivative is the effective electronegativity of atom \(i\), \(-\frac{\partial U}{\partial q_i} = \chi_i\). Thus, the effective force driving the change in atomic charges is the electronegativity difference from the avarage
In the damped dynamic equilibration, the charges are developed dynamically according to the equation of motion with an added damping (friction) term \(- \eta \dot{q}_i\)
This leads to the charges being driven towards a state where the driving force vanishes \(\Delta \chi_i = 0\), i.e., the electronegativities are equal.
During simulation such as molecular dynamics or geometry optimization, charge equilibration is done by running the damped charge dynamics (1) before each force or energy evaluation.
Keywords:
>>> names_of_parameters('dynamic')
['n_steps', 'timestep', 'inertia', 'friction', 'tolerance']
Similarly to Damped dynamics, a potentiostat drives the atomic charges \(q_i\) with the dynamics
Charge equilibration tries to optimize the charge distribution (with fixed total charge) so that the total energy is minimized. The potentiostat, however, connects the system to an electrode (a charge resevoir) at constant potential \(\Phi\). This means that the total charge of the system is allowed to change so that the electronegativity of each atom equals the external potential \(\chi_i = -\frac{\partial U}{\partial q_i} = \Phi\).
Keywords:
>>> names_of_parameters('potentiostat')
['n_steps', 'timestep', 'inertia', 'potential']
This charge optimization method uses the Scipy fmin_slsqp function for constrained optimization using Sequential Least Squares Programming (SLSQP).
Instead of dynamic simulation, the algorithm searches for energy minimum \(U(\mathbf{q})\) with the constraint of constant charge \(\sum_i q_i = Q\).
As the external function does not support callback functions, observers cannot be attached to this algorthm.
Keywords:
>>> names_of_parameters('optimize')
['n_steps', 'tolerance']
Below is a list of methods in ChargeRelaxation, grouped according to the type of functionality.
A class for handling charge dynamics and relaxation.
Pysic does not implement molecular dynamics or geometric optimization since they are handled by ASE. Conceptually, the structural dynamics of the system are properties of the atomic geometry and so it makes sense that they are handled by ASE, which defines the atomic structure in the first place, in the ASE Atoms class.
On the other hand, charge dynamics are related to the electronic structure of the system. Since ASE is meant to use methods such as density functional theory (DFT) in the calculators is employs, all electronic properties are left at the responsibilty of the calculator. This makes sense since in DFT the electron density is needed for calculations of forces and energies.
Pysic is not a DFT calculator and there is no electron density but the atomic charges can be made to develop dynamically. The ChargeRelaxation class handles these dynamics.
Parameters:
Attach a callback function.
Call the given function with the given arguments at constant intervals, after a given number of relaxation steps.
Parameters:
Calls all the attached callback functions.
Parameters:
Performs the charge relaxation.
The relaxation is always performed on the system associated with the Pysic calculator joint with this ChargeRelaxation. The calculated equilibrium charges are returned as a numeric array.
If an ASE Atoms structure is known by the ChargeRelaxation (given through set_atoms()), the charges of the structure are updated according to the calculation result. If the structure is not known, the charges are updated in the structure stored in the Pysic calculator, but not in any other object. Since Pysic only stores a copy of the structure it is given, the original ASE Atoms object will not be updated.
Removes all observers.
Equivalent to set_relaxation_pipe(None).
Returns the atoms object known by the algorithm.
This is the ASE Atoms which will be automatically updated when charge relaxation is invoked.
Returns the Pysic calculator assigned to this ChargeRelaxation.
Return the piped ChargeRelaxation objects.
By default, a list is given containing all the ChargeRelaxation objects in the pipeline (or an empty list).
If called with the argument False, only the ChargeRelaxation next-in-line will be returned (or None).
Parameters:
Creates a dictionary of parameters and initializes all values to 0.0.
Names of the charge relaxation algorithms available.
These are keywords needed when creating the ChargeRelaxation objects as type specifiers.
Short descriptions of the relaxation parameters.
Names of the parameters of the charge relaxation algorithms.
Lets the relaxation algorithm know the atomic structure to be updated.
The relaxation algorithm always works with the structure stored in the Pysic calculator it knows. If pass_to_calculator = True, this method also updates the structure known by the calculator. However, this is not the main purpose of letting the ChargeRelaxation know the structure - it is not even necessary that the structure known by the relaxation algorithm is the same as that known by the calculator.
The structure given to the algorithm is the structure whose charges it automatically updates after relaxing the charges in charge_relaxation(). In other words, if no structure is given, the relaxation will update the charges in the strucure known by Pysic, but this is always just a copy and so the original structure is left untouched.
Parameters:
Assigns a Pysic calculator.
The calculator is necessary for calculation of electronegativities. It is also possible to automatically assign the charge relaxation method to the calculator by setting reciprocal = True.
Note though that it does make a difference whether the calculator knows the charge relaxation or not: If the Pysic has a connection to the ChargeRelaxation, every time force or energy calculations are requested the charges are first relaxed by automatically invoking charge_relaxation(). If there is no link, it is up to the user to start the relaxation.
Parameters:
Sets a given parameter to the desired value.
Parameters:
Sets the numeric values for all parameters.
Use get_parameters() to see the correct order of values.
Parameters:
Sets the numeric values for all parameters.
Equivalent to set_parameter_values()
Parameters:
Sets the relaxation method.
The method also creates a dictionary of parameters initialized to 0.0 by invoking initialize_parameters().
Parameters:
Add another ChargeRelaxation to take place immediately after this relaxation.
Sometimes it is most efficient to start the relaxation with a robust but possibly slow algorithm and switch to a more advanced one once closer to equilibrium. On the other hand, when relaxation is needed between MD steps, this is not possible manually. Therefore, one can link several ChargeRelaxation objects together to automatically invoke a sequential relaxation procedure.
You should only give one ChargeRelaxation for piping. If you want to pipe more than two relaxation sequences together, you should always link the additional relaxation to the currently final ChargeRelaxation object to create a nested pipeline:
rel1 = ChargeRelaxation()
rel2 = ChargeRelaxation()
rel3 = ChargeRelaxation()
rel1.set_relaxation_pipe(rel2)
rel2.set_relaxation_pipe(rel3)
rel1.charge_relaxation() # will do rel1 -> rel2 -> rel3
Parameters: